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1  Executive  Summary 

Multi-material  interfaces  are  sites  of  failure  initiation  in  composite  materials.  The  develop¬ 
ment  of  reliable  quantitative  criteria  for  failure  initiation  in  electronic  components,  adhesively 
bonded  joints  and  laminated  composites  is  obviously  very  important.  At  present  there  are  no 
universally  accepted  procedures  for  the  evaluation  of  fatigue  and  durability  characteristics 
of  structural,  mechanical  and  electronic  components  made  of  composite  materials. 

This  project  was  concerned  with  the  development  of  mathematical  methods  and  compu¬ 
tational  procedures  for  the  determination  of  functionals  which  can  be  correlated  with  failure 
initiation  events  in  composite  materials  subjected  to  thermal  and  mechanical  loads.  The 
approach  is  based  on  the  assumption  that  failure  initiation  events  are  associated  with  the 
natural  straining  modes,  analogously  to  the  well  established  correlation  between  generalized 
stress  intensity  factors  in  linear  elastic  fracture  mechanics  and  crack  propagation  events. 

Failure  criteria  must  be  formulated  in  terms  of  functionals  the  exact  values  of  which  are 
finite.  Stresses  corresponding  to  the  exact  solution  are  usually  infinity  in  singular  points.  The 
numerically  computed  stresses  are,  of  course,  finite  but  very  sensitive  to  the  discretization. 
Therefore  stresses  cannot  be  used  for  formulating  failure  criteria. 

The  specific  objectives  of  the  project  were: 
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1.  Develop  procedures  for  the  numerical  determination  of  the  eigenpairs  A*  and  <f>i  that 
characterize  the  natural  straining  modes  and  natural  flux  states  at  singular  points  in 
heterogeneous  bodies. 

2.  Develop  a  method  for  numerical  determination  of  the  generalized  stress  intensity  and 
flux  intensity  factors  in  heterogeneous  bodies  subjected  to  thermal  and  mechanical 
loading. 

In  addition,  development  of  methods  for  the  estimation  of  limit  loads  for  fiber-reinforced 
composites  in  compression  was  undertaken.  The  instability  of  fibers  is  an  important  consid¬ 
eration  in  the  design  of  fiber-reinforced  composites. 


1.1  Summary  of  accomplishments 

•  A  reliable  numerical  method  for  the  determination  of  the  flux  and  stress  fields  at 
multi-material  interfaces  in  thermoelastic  problems  has  been  developed.  The  method 
involves  numerical  determination  of  the  eigenpairs  of  the  asymptotic  expansion  by  a 
procedure  called  the  modified  Steklov  method  and  determination  of  the  coefficients  by 
a  procedure  based  on  the  complementary  energy  principle. 

•  A  test  implementation  has  been  completed  and  the  effectiveness  of  the  method  estab¬ 
lished  through  benchmark  studies. 

•  Industrial  application  has  been  made  possible  by  a  Phase  I  STTR  grant  to  Engineering 
Software  Research  and  Development,  Inc.  (ESRD)  and  Washington  University.  The 
two-dimensional  thermoelastic  capabilities  have  been  implemented  in  the  commercial 
FEA  code  StressCheck.  This  is  an  unique  capability  which  is  now  being  used  for 
investigation  of  the  correlation  of  observed  failure  events  in  lap-shear  test  specimen 
with  generalized  stress  intensity  factors.  This  work,  started  on  August  1,  1997  is  being 
performed  in  collaboration  with  Raytheon  TI  Systems. 

•  The  original  scope  of  work  was  extended  to  include  numerical  simulation  of  failure  of 
homogeneous  and  composite  elastic  materials  through  loss  of  stability.  This  work,  per¬ 
formed  in  collaboration  with  Dr.  Manil  Suri  of  the  University  of  Maryland  and  Dr.  Ivo 
Babuska  of  The  University  of  Texas,  Austin,  led  to  a  clarification  of  some  fundamen¬ 
tal  theorems  related  to  the  numerical  simulation  of  problems  of  elastic  stability.  A 
doctoral  dissertation  has  been  completed  [3]. 
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1.2  Personnel  supported 

Faculty: 

Dr.  Barna  A.  Szabo  (PI) 

Post-doctoral  persons: 

Dr.  Ricardo  L.  Actis  (part  time) 

Dr.  Xian-Zhong  Guo  (part  time) 

Dr.  Gyorgy  Kiralyfalvi 
Graduate  Students: 

Mr.  Andre  Tamagnini  Noel  (D.Sc.  candidate,  Graduated  May,  1996) 
Mr.  Gyorgy  Kiralyfalvi  (D.Sc.  candidate,  Graduated  May,  1997) 

Ms.  Li  Zhang  (D.Sc.  candidate) 


1.3  Consultative  and  advisory  functions 

The  Principal  Investigator  presented  briefings  to  eight  Air  Force  contractors  and  one  Navy 
laboratory: 

1.  McDonnell  Douglas  Aerospace  St.  Louis,  MO  (contact  persons:  Mr.  Scott  Fields,  Mr. 
Daniel  Dudley)  April  25,  1996 

2.  Boeing  Aerospace,  Downey,  CA.  Dr.  Saeed  Paydarfar  visited  Washington  University 
and  was  briefed  on  the  scope  and  objectives  of  the  project.  May  19,  1997 

3.  Boeing  Aircraft  Co.,  Wichita,  KS  (contact  person:  Mr.  Phillip  Legate)  June  4,  1997 

4.  Cessna  Aircraft  Co,  Wichita,  KS  (contact  person:  Mr.  Milan  Radovanov)  June  5,  1997 

5.  Raytheon  TI  Systems,  Dallas,  TX  (contact  person:  Dr.  Terry  Baughn)  July  29,  1997 

6.  Lockheed-Martin,  Fort  Worth,  TX  (contact  person:  Mr.  Michael  Barnhart)  July  30, 
1997 

7.  Allied  Signal,  Phoenix,  AZ  (contact  person:  Dr.  Malak  Malak  602-231-3701)  Septem¬ 
ber  5,  1997 

8.  Structures  Division,  NAWCAD,  Patuxent  River,  MD  (contact  person:  Dr.  David  Bar¬ 
rett  301-342-9360)  September  8,  1997 
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9.  Lockheed-Martin,  Marietta,  GA  (contact  person:  Dr.  Stephen  P.  Engelstad  770- 
494-9714)  March  3,  1998.  This  presentation  was  made  to  members  of  a  govern¬ 
ment/industry  consortium  known  as  the  Composites  Affordability  Initiative.  Members 
represent  each  of  the  major  US  aerospace  companies  and  the  Air  Force  and  Navy.  This 
particular  meeting  was  hosted  by  Lockheed  Martin. 


1.4  Publications  and  presentations 


[1]  Noel,  A.  and  B.  Szabo.  Formulation  of  geometrically  non-linear  problems  in  the  spatial 
reference  frame.  Int.  J.  Numer.  Methods  Eng .,  40:1263-1280,  1997. 

[2]  Noel  A.  T.  Spatial  Formulation  and  Numerical  Solution  of  Geometrically  Nonlinear 
Problems  in  Finite  Elasticity.  D.Sc.  Dissertation,  Sever  Institute  of  Technology,  Wash¬ 
ington  University,  St.  Louis,  Missouri,  1996. 

[3]  I.  Babuska,  and  B.  Szabo.  New  problems  and  trends  in  the  finite  element  method.  In 
J.  R.  Whiteman,  editor,  The  Mathematics  of  Finite  Elements  and  Applications ,  pages 
1-33,  Chichester,  1997.  John  Wiley  and  Sons. 

[4]  B.  Bertoti,  E.  and  Szabo.  Adaptive  selection  of  polynomial  degrees  on  a  finite  element 
mesh.  To  appear  in  Int.  J.  Numer.  Meth.  Engng.,  1998. 

[5]  I.  Paczelt  and  T.  Szabo,  B.  and  Szabo.  Solution  of  elastic  contact  problems  by  the 
p-version  of  the  finite  element  method.  4th  U.S.  National  Congress  on  Computational 
Mechanics,  August,  1997. 

[6]  S.  A.  Prost-Domaski,  B.  A.  Szabo,  and  G.  I.  Zahalak.  Large-deformation  analysis  of 
nonlinear  elastic  fluids.  Computers  and  Structures ,  64:1281-1290,  1997. 

[7]  B.  Szabo  and  Z.  Yosibash.  Numerical  analysis  of  singularities  in  two  dimensions.  Part  2: 
Computation  of  generalized  flux/stress  intensity  factors.  Int.  J.  Numer.  Meth.  Engng., 
39:409-434,  1996. 

[8]  B.  Szabo  and  Z.  Yosibash.  Superconvergent  extraction  of  flux  intensity  factors  and  first 
derivatives  from  finite  element  solutions.  Comput.  Meth.  Appl.  Mech.  Engrg.,  129:349- 
370,  1996. 

[9]  B.  A.  Szabo.  Hierarchic  models  and  discretizations.  Symposium  on  Advances  in  Com¬ 
putational  Mechanics,  The  University  of  Texas  at  Austin,  January  1997. 
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[10]  B.  A.  Szabo  and  R.  A.  Actis.  Failure  analysis  of  composite  materials.  International 
Mechanical  Engineering  Congress  and  Exposition,  November,  1996  Atlanta,  GA. 

[11]  G.  Kiralyfalvi.  Linear  Models  of  Elastic  Stability.  D.Sc.  Dissertation,  Sever  Institute  of 
Technology,  Washington  University,  St.  Louis,  Missouri,  1997. 

[12]  B.  A.  Szabo  and  G.  Kiralyfalvi.  Mathematical  models  of  buckling  and  stress  stiffening. 
4th  U.S.  National  Congress  on  Computational  Mechanics,  August,  1997. 

[13]  R.  Szabo,  B.  and  Actis.  The  problem  of  model  selection  in  numerical  simulation.  4th 
U.S.  National  Congress  on  Computational  Mechanics,  August,  1997. 

[14]  Y.  Volpert,  T.  Szabo,  I.  Paczelt,  and  Szabo  B.  Application  of  the  space  enrichment 
method  to  problems  of  mechanical  contact.  Finite  Elements  in  Analysis  and  Design, 
24:157-170,  1997. 

[15]  Y.  Wang,  P.  Monk,  and  B.  Szabo.  Computing  cavity  modes  using  the  p-version  of  the 
finite  element  method.  IEEE  Transactions  of  Magnetics,  32:1934-1940,  1996. 

[16]  Z.  Yosibash  and  B.  A.  Szabo.  A  note  on  numerically  computed  eigenfunctions  and 
generalized  stress  intensity  factors  associated  with  singular  points.  Engineering  Fracture 
Mechanics,  54(4):593-595,  1996. 

[17]  Z.  Yosibash  and  B.  A.  Szabo.  Failure  analysis  of  composite  materials  and  multi-material 
interfaces  Proceedings,  1995  Design  Engineering  Technical  Conferences,  ASME  DE- 
Vol.  83:133-139,  1995. 

[18]  Z.  Yosibash  and  B.  Szabo.  Numerical  analysis  of  singularities  in  two  dimensions.  Part 
1:  Computation  of  eigenpairs.  Int.  J.  Numer.  Meth.  Engng.,  38:2055-2082,  1995. 

[19]  Z.  Yosibash  Numerical  thermo-elastic  analysis  of  singularities  in  two-dimensions.  In¬ 
ternational  Journal  of  Fracture,  74:341-361,  1996. 

[20]  Z.  Yosibash  Computing  edge  singularities  in  elastic  anisotropic  three-dimensional  do¬ 
mains  International  Journal  of  Fracture,  86:221-245,  1997. 

[21]  G.  Kiralyfalvi  and  B.  Szabo.  Quasi-Regional  Mapping  for  the  p-Version  of  the  Finite 
Element  Method.  Finite  Elements  in  Analysis  and  Design ,  27:85-97  1997. 


1.5  Transitions 

The  new  capabilities  have  been  made  available  to  Air  Force  laboratories  and  contractors 
through  a  professional  quality  software  called  StressCheck.  StressCheck  is  being  developed 
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and  marketed  by  Engineering  Software  Research  and  Development,  Inc.,  located  in  St.  Louis, 
MO.  The  current  users  of  Stress  Check  include  Boeing  Aircraft  Company  (on  several  loca¬ 
tions);  Piper  Aircraft  Co.,  Northrop  Grumman  Corporation;  Cessna  Aircraft  Co.  NASA 
Johnson  Space  Center  and  others. 

A  collaborative  effort  was  started  with  Raytheon  TI  System  for  an  experimental  inves¬ 
tigation  of  realtionships  between  generalized  stress  intensity  factors  and  failure  initiation 
events  at  bonded  interfaces. 

McDonnell  Douglas  Aerospace  St.  Louis  (now  Boeing)  funded  a  project  with  ESRD  for  a 
particular  specialization  of  the  material  and  geometric  nonlinear  analysis  capabilites  within 
the  p-version  of  the  finite  element  method  for  application  to  the  analysis  of  cold-worked 
holes  and  attachment  lugs.  This  technology  was  developed  at  Washington  University  under 
AFOSR  sponsorship. 

2  Technical  description 

Singular  points  are  those  points  in  a  structural  component  where  a  reentrant  corner  occurs 
(like  cracks  and  V-notches),  material  properties  abruptly  change  along  a  free  edge,  interior 
points  of  three  (or  more)  zones  of  different  materials  intersect,  or  an  abrupt  change  in 
boundary  conditions  occurs,  see  Fig.  1. 


Interior  Reentrant  Edge 


Figure  1:  Typical  singular  points  associated  with  multi-material  interfaces. 

In  the  vicinity  of  these  singular  points  the  exact  solution  of  the  problem  of  elasticity  is 
of  the  form: 

OO 

uex  =  Airx,<j>i(6)  (1) 

t=i 

where  uex  is  the  exact  displacement  vector  function,  r  and  0  are  polar  coordinates  centered 
on  the  singular  point,  4>i  =f  {<f>i\x  4>i\y}  is  a  piecewise  smooth  vector  function,  and  A;  are 
coefficients.  Eq.  (1)  is  an  asymptotic  expansion  of  the  exact  solution  at  the  singular  point. 
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Within  a  radius  of  convergence  the  exact  solution  of  the  problem  of  elasticity  can  be  written 
in  this  form.  The  exponents  A;  (numbered  such  that  Ai  <  A2  <  A3 . . .)  and  the  corresponding 
functions  fa(0)  depend  on  the  material  properties  and  the  geometric  details  at  the  singular 
point.  These  can  be  determined  by  solving  an  eigenvalue  problem.  Details  are  given  in 
Section  2.3.  A  well  known  example  is  linear  elastic  fracture  mechanics  in  two  dimensions 
where  A:  =  1/2  and 


H 

II 

1 

17 

1  30 

2 G 

[v  ~  2y 

(cos-- 

-  cos  — 
2  2 

^1|  y 

1 

[/  1> 

l  .  0 

1  .  3 O' 

2  G 

Vs  +  2; 

(sin-- 

—  sm  — 

2  2 
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G  is  the  shear  modulus,  k  =  3  —  4u  for  plane  strain,  k  =  (3  —  v)/(l  +  v)  for  plane  stress 
where  u  is  Poisson’s  ratio. 


In  linear  elastic  fracture  mechanics  A,-  and  fa  have  been  determined  by  classical  methods, 
and  only  the  stress  intensity  factor,  which  is  proportional  to  Ai,  has  to  be  determined  by 
numerical  means.  For  details  see,  for  example,  [6]. 

In  the  general  case  of  multi-material  singularities,  such  as  those  shown  in  Fig.  1,  not 
only  A(  but  also  A,-  and  fa{6)  have  to  be  determined  by  numerical  means.  The  elastic  stress 
is  infinity  in  the  singular  point  when  0  <  Ai  <  1  and  Ai  ^  0  and/or  0  <  A2  <  1  and 
A2  7^  0,  etc.  A  natural  straining  mode  is  the  strain  state  associated  with  a  particular 
term  of  the  asymptotic  expansion,  eq.  (1).  As  explained  in  Section  2.1,  the  natural  straining 
modes  provide  a  linkage  between  linear  computations  and  observed  failure  initiation  or  failure 
propagation  events. 

In  heat  conduction  the  asymptotic  expansion  is  analogous  to  eq.  (1): 

OO 

rH  =  EV#)  (4) 

i=l 

where  Tex  is  the  exact  solution  of  the  heat  conduction  problem  and  fa  are  piecewise  smooth 
scalar  functions.  The  coefficients  A,  are  called  flux  intensity  factors. 

Numerical  accuracy  is  essential  because  unless  the  accuracy  of  the  computed  data  is 
known  it  would  not  be  possible  to  tell  whether  the  working  hypothesis  is  wrong  or  the 
numerical  errors  are  too  large,  or  both.  In  some  cases  a  large  error  in  the  working  hypothesis 
is  nearly  canceled  by  a  similarly  large  numerical  error,  leading  to  false  conclusions. 


2.1  Basic  principles  and  assumptions 

Consider  the  neighborhood  of  a  singular  point.  It  is  assumed  in  the  following  that  the 
principles  of  continuum  mechanics  remain  valid  everywhere  within  the  body  up  to  the  fail- 
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ure  initiation  event.  The  possibility  of  strongly  nonlinear  behavior  in  the  neighborhood  of 
singular  points  is  not  excluded,  however. 


Figure  2:  Definition  of  Tnl  and  Tel- 

Let  unl  =  {ux  uy}NL  be  the  solution  of  the  general  nonlinear  continuum  mechanics 
problem.  It  is  expected  that  failure  initiation  will  depend  on  unl  (more  precisely  some 
functionals  computable  from  unl)  in  the  strongly  nonlinear  region  of  the  singular  point 
bounded  by  a  boundary  Tnl,  as  shown  in  Fig.  2.  This  region  is  called  the  process  zone.  Let 
Tel  be  a  curve  outside  of  Tnl,  and  let  G  be  an  operator  which  associates  the  solution  unl 
of  the  nonlinear  problem  inside  r./vx  with  the  boundary  condition  qel  specified  on  Tel,  that 
is: 

G(gEL )  —  uex  ,  9el  =  UNL\rEL 

where  §el  —  unl\vel  denotes  the  trace  of  Hex  on  Tel-  Denote  the  exact  solution  of  the 
linear  elastic  problem  by  uel -  The  basic  assumptions  (which  are  valid  in  linear  elastic 
fracture  mechanics)  are  stated  in  the  following: 

Assumption  A: 

Inside  of  Tel  the  error  G(unl\tbl)  ~  G{uel\tel)  is  so  small  that  conclusions  based  on 
G(uel\tel)  are  sufficiently  close  to  the  conclusions  based  on  G{u^e\tel)  for  practical  pur¬ 
poses.  This  assumption  is  expected  to  be  valid  whenever  the  nonlinear  behavior  is  confined 
entirely  to  some  small  region  inside  Tel- 

Assumption  A  leads  to  the  important  conclusion  that  failure  initiation,  which  depends  on 
the  solution  of  the  nonlinear  problem  inside  of  T^l,  can  be  determined  through  the  solution 
of  the  linear  elastic  problem,  even  though  all  basic  assumptions  of  the  linear  theory  may  be 
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violated  inside  of  T/vl-  Consequently  failure  initiation  in  the  neighborhood  of  the  singular 
point  can  be  predicted  on  the  basis  of  the  linear  theory  of  elasticity.  (This  is  because  uel 
defines  uNL\rBL-) 

Assumption  B: 

There  exists  a  physical  principle  which  establishes  the  relationship  between  crack  initi¬ 
ation  and  the  stress  field  on  the  basis  of  information  obtained  from  the  linear  solution  uel 
only.  Linear  elastic  fracture  mechanics  (LEFM)  is  a  typical  application  of  Assumption  B. 

In  general,  the  linear  solution  uel  is  not  known,  only  an  approximation  to  uel,  which 
will  be  denoted  by  ufe,  is  known.  Therefore  the  following  assumption  is  necessary: 

Assumption  C: 

There  exist  a  norm  |j  •  ||  such  that  when  \\uel  —  ufe\\  is  sufficiently  small  then  the  physical 
principle  of  Assumption  B  is  not  sensitive  to  the  replacement  of  uel  with  ufe-  Of  course, 
the  norm  ||  •  ||  is  expected  to  depend  on  the  physical  principle  of  Assumption  B,  which  is 
material-dependent.  If  conclusions  are  to  be  based  on  ufe  then  ufe  has  to  be  close  to  uel 
in  this  particular  norm. 

Based  on  these  assumptions  linear  computations  can  be  used  for  the  prediction  of  fail¬ 
ure  initiation  and  failure  propagation  even  though  failure  processes  are  highly  nonlinear  in 
nature.  There  are  two  key  elements  of  failure  initiation  analysis: 

1.  A  hypothesis  concerning  the  relationship  between  certain  parameters  of  the  stress  or 

strain  field  and  observed  failure  initiation  or  crack  propagation  events. 

2.  Convincing  experimental  confirmation  that  the  hypothesis  holds  independently  of  vari¬ 

ations  in  geometric  attributes,  loading  and  constraints. 

It  would  not  sensible  to  perform  experiments  without  a  hypothesis  based  on  the  functionals 
that  characterize  the  stress  or  strain  fields  in  the  neighborhood  of  critical  points  and  com¬ 
putations  cannot  provide  useful  information  about  the  conditions  under  which  failure  occurs 
without  experimental  data. 

For  details  on  the  algorithms  developed  for  the  computation  of  the  natural  straining 
modes  at  multi-material  interfaces  and  the  generalized  stress  intensity  factors  we  refer  to 
[12],  [7],  [8]  [10],  [11]. 

Remark  2.1  The  assumption  that  the  material  is  elastic  on  Tel  is  not  essential.  Similar 
considerations  apply  to  nonlinear  elasticity  and  the  deformation  theory  of  plasticity.  In  fact, 
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the  methods  of  LEFM  have  been  extended  to  the  deformation  theory  of  plasticity  through 
the  use  of  the  J-integral  [2],  [5]. 

In  the  following  the  procedures  developed  for  the  computation  of  eigenpairs  and  their 
coefficients  in  heat  conduction  and  elasticity  are  outlined  and  illustrated  by  examples.  Ad¬ 
ditional  information  can  be  obtained  from  the  references  listed. 


2.2  The  problem  of  heat  conduction 

The  index  notation  is  used  in  the  following.  For  two  dimensional  problems  the  range  of  the 
indices  is  2,  and  for  three  dimensional  problems  the  range  is  3.  The  summation  convention  is 
used.  The  formulation  of  the  mathematical  is  described  for  the  problem  of  heat  conduction. 

The  heat  balance  equation  is  analogous  to  the  equilibrium  equation  in  elasticity: 

-qi,i  +  <2  =  0  (5) 

where:  qi  is  the  flux  vector  (in  W/m2  units)  and  Q  is  the  rate  of  heat  generation  per  unit 
volume  (in  W/m3  units).  Multiplying  eq.  (5)  by  the  scalar  function  v,  integrating  and 
applying  Green’s  lemma,  we  have: 


Using  Fourier’s  law  of  heat  conduction  (which  is  analogous  to  Hooke’s  law): 

qi  =  Ktj'T'j 

where  Kij  is  assumed  to  be  independent  of  T,  we  have  for  all  v  €  Hl(Q): 

f  KijVtiT,jdV  =  /  QvdV-  [  qimvdS.  (6) 

JQ,  vF 

This  is  analogous  to  the  principle  of  virtual  work  in  elasticity.  Alternatively,  the  exact 
solution  of  the  heat  conduction  problem  is  the  minimizer  of  the  functional  tt: 

t r(T)  =  I-  f  KitfiTj  dV  -  f  QTdV  +  I  qimTdS  (7) 

Li  vii  j o  jt 

on  the  set  of  the  admissible  temperature  fields.  This  is  analogous  to  the  principle  of  minimum 
potential  energy  in  elasticity. 
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Figure  3:  Typical  singular  point  in  2D.  Definition  of  Hr. 


2.3  Computation  of  A*  and  4>i(9) 


The  procedure  for  the  numerical  approximation  of  the  eigenpairs  A,  and  <j>i(0)  is  briefly 
described  for  the  problem  of  heat  conduction.  Additional  details  are  available  in  [8]. 

Consider  a  small  neighborhood  of  a  singular  point  bounded  by  a  circle  of  radius  R  and 
denoted  by  Hr,  as  shown  in  Fig.  3.  Let  Q  =  0  on  Hr  x  t  where  t  is  the  thickness,  assumed 
to  be  constant,  and  «  =  0  on  Ti  and  T2  and  seek  solutions  of  the  form: 

T  =  rx<f>{9). 


Noting  that 

dr  dr  dr  do  dr  Q  1  dr  .  a 

a  a  +  Q  =  —  COS  0 - TrSintl 

dr  dx i  dd  dxi  dr  r  dO 

dr  dr  dr  do  dr  .  a  idr  Q 

a - 1"  =  T-sinH  -^r-cosO 

dr  dx2  d6  dx 2  dr  r  dO 

and  nl  =  {cos  6  sin  0],  we  write: 


We  are  now  in  a  position  to  apply  eq.  (6)  which  results  in: 


{B(r,  V )  -  A /XT,  v))  -  A M(r,  v)  =  0 


(8) 


where: 


B(7»  [  Kijv/TjiV 

M(T,v )  =  /  (Kn  cos2  9  +  K12  sin  29  +  K22  sin2  9)T v  td.0 

Jtr 


M{T,  v)  =f  f  ((K22  — Ku)  sin 0 cos $  + Ku  cos 29) -^jrvtdO 
JrR  09 

We  seek  A  >  0,  T  €  R^Dr),  dT/00  G  L2(Tr)  such  that  eq.  (8)  is  satisfied  for  all  v  e 
This  non-symmetric  eigenvalue  problem  can  be  solved  numerically. 

Remark  2.2  In  the  case  of  isotropic  materials  N{T,  v)  =  0  hence  the  eigenvalue  problem 
is  symmetric.  In  the  special  case  Kij  =  KSij ,  where  K  is  constant,  we  have: 


/  TiVidV-Xf  TvtdO  =  0. 


The  corresponding  strong  form  is: 
subject  to  the  boundary  conditions 


AT  =  0 


OT  X 

u  =  0  on  Ti,  r2;  — T  on  Tr. 

Or  R 

Remark  2.3  Af(T,  v)  is  non-symmetric,  nevertheless  all  eigenvalues  are  real. 

Remark  2.4  A,-  do  not  converge  monotonically.  No  minimum  principle  is  involved. 
Remark  2.5  The  formulation  is  analogous  for  any  set  of  homogeneous  boundary  conditions 

on  rx,  r2. 


The  Steklov  method  on  Qr  requires  hp-meshing.  This  is  because  the  rate  of  convergence 
of  p-extensions  is  low  due  to  the  presence  of  the  singular  point.  Therefore  it  is  better  to 
use  a  modified  domain  f shown  in  Fig.  4.  Using  Q,r  is  called  the  modified  Steklov  method. 
Detailed  discussions  on  the  procedures  are  available  in  references  [8],  [12],  [7]. 

Remark  2.6  The  modified  Steklov  method  will  find  not  only  those  eigenfunctions  which  lie 
in  /T(fttf)  but  also  eigenfunctions  corresponding  to  negative  values  of  A. 
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Figure  4:  The  domain  fljEj. 


2.4  Extraction  of  the  flux  intensity  factors 

Analogously  to  eq.  (1)  the  temperature  field  in  the  neighborhood  of  singular  points  is  of  the 
form: 

OO 

T=Y,Ay~<t>m  (9) 

%- 1 

where:  A,-  represents  the  generalized  £ux  intensity  factors,  A;  and  4>i  are  the  eigenpairs 
characterized  by  the  topological  details  at  the  singular  point  and  the  material  properties. 

Once  the  finite  element  solution  is  available,  the  flux  intensity  factors  can  be  computed 
from  the  finite  element  solution  by  the  contour  integral  method,  the  complementary  energy 
method,  and  the  T2  projection  method.  These  are  briefly  described  in  the  following. 

2.4.1  The  contour  integral  method 

The  contour  integral  method  is  a  procedure  which  utilizes  the  orthogonality  of  eigenfunctions 
and  the  fact  that  if  A*  is  an  eigenfunction  then  —A;  is  also  an  eigenfunction.  For  details  we 
refer  to  [1],  [6]. 

2.4.2  The  complementary  energy  method 

Define: 

nc(ft)=fo/  CiiWdV-  f  qimTFEtde 

Ci  JQji  «/r ji 

where  C,j  is  the  inverse  of  Kif,  q, ;  satisfies  the  heat  balance  equation:  —qi,i  +  (5  =  0  (and 
prescribed  homogeneous  flux  boundary  conditions);  Tfe  is  the  temperature  computed  by 
the  finite  element  method. 


Letting  q,-  =  —KijT,j  and  minimizing  the  complementary  energy  functional  IIc(q,)  with 
respect  to  A,,  yields  approximations  to  the  generalized  flux  intensity  factors. 

2.4.3  The  L2  projection  method 

This  method  involves  the  projection  of  ufe  onto  the  space  spanned  by  the  eigenfunctions. 
Specifically,  A,  are  determined  from  the  condition: 

jf  -  X]  AirXi <})i(9) j  rXj(j>j(0)  rdrdO  =  0  j  —  1, 2, . . . ,  n 

which  yields  n  equations  for  Ai,  i  =  1, 2, . . . ,  n. 


2.5  Example:  The  slit  domain  problem 

Consider  the  problem  AT  =  0  on  a  unit  circle  slit  along  the  axis  x2  with  the  boundary 
conditions 

T  =  0  on  Ti,  q2  =  0  on  T2,  qn  =  q;n»  =  x2  on  r3. 

In  this  case:  Ai  =  0.25,  A2  =  0.75,  A3  =  1.25.  This  is  a  very  challenging  problem  from  the 
point  of  view  of  numerical  approximation  by  the  finite  element  method,  owing  to  the  fact 
that  the  lowest  eigenvalue  is  1/4,  hence  the  theoretical  rate  of  convergence1  of  the  p- version 
is  1/4  and  the  theoretical  rate  of  convergence  of  the  h- version  is  1/8.  The  error  is  most 
effectively  controlled  by  hp-extension,  utilizing  geometrically  graded  meshes  [6]. 

The  finite  element  mesh,  consisting  of  12  elements,  and  the  temperature  distribution  at 
p  =  8  (trunk  space)  is  shown  in  Fig.  6.  The  p-convergence  of  the  energy  and  A\,  A2,  As  are 
shown  in  Table  1.  The  extrapolated  values  are  shown  in  the  last  row.  It  is  seen  that  the 
rate  of  convergence  is  close  to  the  theoretical  rate  of  0.25. 


2.6  Example:  Two-material  internal  interface  problem 

Consider  the  problem  AT  =  0  on  a  unit  circle.  On  the  first  quadrant  has  the  material 
properties  are  Kn  =  K22  =  10.0,  K\2  =  0  on  the  other  three  quadrants  Ku  =  K22  =  1.0, 
K12  =  0.  On  the  boundary  qn  =  f(6)  where  f(0)  is  a  function  given  by  Oh  and  Babuska  in 

[4]- 

xThe  theoretical  rate  of  convergence  is  given  by  /?  in  the  a  priori  estimate  \\uex  —  «F£|U(n)  <  kN~f3 
were  k  is  a  positive  constant  and  N  is  the  number  of  degrees  of  freedom.  See  [6]. 
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Figure  5:  The  slit  domain  problem. 


2.7  The  problem  of  thermo-elasticity 

The  equations  of  equilibrium  are: 

0* j,j  “h  F1,-  =  0  (10) 

Multiplying  (10)  by  a  test  function  vt  and  applying  Green’s  lemma,  the  generic  form  of  the 
principle  of  virtual  work  is  obtained: 

j  dV  =  J  FiVi  dV  +  jf  TX-  dS 

where 

4?  =  \ivi,3  +  vi,i) 

is  the  small  strain  tensor  corresponding  to  the  virtual  displacement  vl  and 

4“ >  =  -  auT) 

is  the  stress  tensor  corresponding  to  a^i  represents  the  coefficients  of  thermal  expansion 
and  T  is  the  temperature. 

Remark  2.7  The  temperature  field  T  is  continuous  but  does  not  have  to  be  continuous. 
Remark  2.8  Eijki  may  be  a  function  of  T. 

An  alternative  formulation  is  the  principle  of  minimum  potential  energy: 

n(«)  = 

—  f  F{Ui  dV  —  f  TiUi  dS 
Jn  Jr 


ID  s  SOL 
Run  =  8 
Fnc.  =  U 

Max  =  1.8981e+0O 
Min  =  -7.9321e-01 

1.9000e+00 

1.6000e+00 

13000e+00 

l.OOOOe+OO 

7.0000e-0I 

4.0000e-01 

1.0000e-01 

-2.0000e-01 

-5.0000e-01 

-8.0000e-01 


Figure  6:  The  slit  domain  problem:  Mesh  layout  and  contour  plot  of  T,  p=8,  trunk  space. 

The  exact  solution  of  the  problem  of  elasticity  minimizes  the  potential  energy: 

n(u£x)  =  min  n(u)  (11) 

u€E(CL) 

where  II  is  the  potential  energy  and  E(  fl)  represents  the  space  of  admissible  functions.  In 
the  case  of  isotropic  materials  the  Euler  equations  are: 

GW  +  (A  +  G)  W),;  =  -Fi  +  pTti 

where 

P  d=  (3A  +  2 G)a. 

where  A  is  the  Lame  parameter,  G  is  the  shear  modulus  and  a  is  the  coefficient  of  thermal 
expansion. 
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Table  1:  p-Convergence  of  the  energy  and  A2,  A3.  Twelve-element  mesh.  The  extrapo¬ 
lated  values  are  shown  in  the  last  row. 


p 

N 

Potential 

Energy 

Rate  of 
Conv. 

Est.’d 
Rel.  Err. 

Ax 

a2 

^3 

1 

12 

-1.977649 

0.00 

35.62 

0.794380 

-0.901532 

0.373088 

2 

36 

-2.141858 

0.39 

23.31 

1.004779 

-0.953011 

0.457475 

3 

66 

-2.178293 

0.29 

19.56 

1.110997 

-0.964365 

0.455113 

4 

108 

-2.196471 

0.24 

17.39 

1.168927 

-0.968574 

0.453430 

5 

162 

-2.208072 

0.23 

15.85 

1.203440 

-0.969836 

0.452681 

6 

228 

-2.216175 

0.22 

14.68 

1.226336 

-0.970035 

0.452696 

306 

-2.222206 

0.22 

13.74 

1.243025 

-0.970053 

0.452696 

396 

-2.226874 

0.22 

12.97 

1.255925 

-0.970047 

0.452703 

OO 

OO 

-2.264978 

0.25 

0 

1.359910 

-0.970047 

0.452696 

2.8  Computation  of  thermal  stress  intensity  factors 

Outline  of  the  solution  algorithm. 


1.  Assuming  that  the  displacement  field  at  the  singular  point  is  of  the  form 

Ui  =  r^i(e) 

compute  the  eigenpairs  fXj,  j  =  1,2,...  using  the  modified  Steklov  method 

where  the  second  index  on  represents  the  ordinal  number  of  the  eigenfunction. 
This  involves  the  solution  of  a  non-symmetric  eigenvalue  problem  of  the  form: 

Vj)  —  Af(ui,  Vi))  -  /j.M(ui,  Vi)  =  0 

The  eigenvalues  are  usually  complex.  Both  the  real  and  imaginary  parts  must  be 
considered. 

(H) 

2.  Construct  the  homogeneous  part  of  the  statically  admissible  stress  field  a  -  ’  from 

<m«) 

3- 1 

where  Cj  represents  the  generalized  stress  intensity  factors.  Specifically, 

/  TJ\  1 

<?ij  =  - Eijki(uk,i  +  uiik). 
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Figure  7:  Problem  definition:  Two-material  internal  interface 


Figure  8:  Two-material  internal  interface  problem:  The  first  four  eigenfunction  (f>i(0),  cor¬ 
responding  to  Ai  =  0.7317. 

3.  The  particular  solution  <r}j  ’  satisfies 

j,j  —  EijkiakiT,j 

on  Hr  and  the  homogeneous  traction  boundary  conditions  on  Ti,  IV  In  general  it 
is  difficult  to  construct  crjp.  Note,  however,  that  crjp  is  of  the  order  rA  where  A  =f 
Amin  >  0  is  the  smallest  eigenvalue  of  the  thermal  problem.  On  the  other  hand,  '  is 
of  order  where  /i  =f  /xm;n  >  0  is  the  smallest  eigenvalue  of  the  elasticity  problem. 

4.  Construct  the  complementary  energy  functional  on  Hr: 

nc(crtj)  =f  \  [  Cijki(Tij<Tki  dV  -  [  (TijUjuf^t  ds 
z  jQr  Jr r 


=  If  ciikl*VV»Jv+[  clit,tfV'>dv 

Z  JSIr  JQji 

v - - - '  - - S, - ' 

0(R2»)  0(HA+^+1) 

+  \  f  CijU^iPcriP  dV  -  I  (TijTijuf^t  ds 

z  Jq,r  JTr 

v - * - ' 

0(R2*+2) 


This  indicates  that  if  R  is  sufficiently  small  then  afjj  may  be  neglected. 

5.  Compute  the  thermal  stress  intensity  factors  C by  minimizing  nc(<r|j?^)  on  Tflfc  for 
a  sequence  of  decreasing  radii  Rk,  k  =  1, 2, . . . , n. 

6.  Use  Richardson  extrapolation  to  find 

Ct  =  lim  Cf\ 

3  Rk->  0  3 

2.9  Example:  Cracked  panel  subject  to  thermal  load 

A  centrally  cracked  panel  is  subjected  to  T  =  100  at  the  perimeter;  T  =  0  on  the  crack 
faces.  Ku  =  K22  =  1-0;  =  0;  E  =  1.0,  u  —  0.3,  a  —  0.01,  plane  strain.  L  =  w  =  50.0, 

a  =  1  is  shown  in  Fig.  9.  The  p-convergence  of  the  first  thermal  stress  intensity  factor  Ai 
computed  with  Richardson  extrapolation,  is  shown  in  Fig.  11  and  the  results  reported  by 
Yosibash  in  [10]  using  direct  computation  are  given  in  Table  2. 


Table  2:  Results  reported  by  Yosibash  using  direct  computation 


R/a 

0.5 

0.1 

0.01 

0.001 

0.0006 

Cx 

0.4908 

0.3781 

0.3528 

0.3491 

0.3481 

2.10  Sources  of  errors 

The  methods  used  for  computing  the  finite  element  solution,  the  eigenpairs  and  the  gener¬ 
alized  flux  intensity  factors,  are  approximate  methods,  hence  certain  errors  are  incurred: 

1.  In  numerical  work  the  asymptotic  expansion  is  truncated  to  a  few  terms. 
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Figure  9:  The  cracked  panel  problem. 

2.  The  eigenpairs  are  only  approximations.  Therefore  qi  do  not  satisfy  the  heat  balance 
equation  exactly. 

3.  The  prescribed  temperature  on  Fr  is  approximate,  computed  from  the  finite  element 
solution. 

Nevertheless,  numerical  experience  has  indicated  that  this  method  of  extraction  is  supercon- 
vergent  well  beyond  the  range  of  accuracy  required  in  engineering  work  [8]. 


2.11  Stability  problems 


This  topic  is  of  substantial  interest  in  aerospace  engineering  because  the  sizes  of  compression 
members  are  determined  primarily  by  stability  considerations.  It  is  also  of  great  importance 
in  micromechanics  where  the  strength  of  the  composite  materials  is  often  determined  by 
the  buckling  of  fibers.  The  formulation  and  investigation  of  stability  problems  in  the  fully 
three-dimensional  setting  was  undertaken.  Details  are  available  in  a  doctoral  dissertation 
[3].  A  brief  outline  is  presented  in  the  following.  Define: 


_o  def  »  * 
aij  ™%j 

where  <r*-  is  the  pre-buckling  stress  state.  We  are  interested  in  finding  u,-  6  E(Cl)  such  that: 

I  CijkiiiijVkj  dV  +  A  /  (7^jUa^va^j  dV  — 

f  FiVi  dV  +  [  Tm  dS  +  /  f  CijkiotkiVij  dV 

J£l  JVj' 
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Figure  10:  The  mesh,  temperature  distribution  and  the  resulting  deformation. 


for  all  Vi  €  E( fl).  Note:  er*-  must  be  such  that: 

|  ( Cijki  T  A<r ji&ik'j'U'ijUkj  dV 


<  oo 


for  all  Ui  €  E(Cl).  The  set  of  A  for  which  a  solution  exists  is  the  resolvent  set.  The  complement 
is  the  spectrum.  The  spectrum  may  be  point,  continuous  or  residual. 

The  work  done  by  the  initial  stress  due  to  the  product  terms  of  the  Green-Lagrange 
strain  tensor  is  incorporated  in  the  strain  energy: 

U(Ui)  =  |  J  Cijkllijtld  dV+\ja  ^a,iUaJ  dV. 

S - V - " 

work  of  o1®  - 

where  Ui  is  a  small  increment  of  displacement  and  e,j  is  the  small  strain. 
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Figure  11:  p-Convergence  of  the  first  thermal  stress  intensity  factor  C\  computed  with 
Richardson  extrapolation. 

Biot  (1938)  and  Prager  (1947)  proposed  a  different  formulation.  Their  definition  is: 

U ( U { )  —  —  f  Cijkl^ij^kl  dV  T  f  €-ai^-aj)  dV 

2  Ja  2  Ja 

The  potential  energy  is: 

II(uj)  =f  U ( Ui )  -  jf  FiUi  dV  —  jf  ^  Titii  dS  +  jf  T' Cijki<Xkiuitj  dV 

We  seek  u;  6  £(fl)  such  that  the  potential  energy  is  stationary: 

XTT/_  ,  def  / (Hi  +  £V{)\  o 

dn(ui)  =  ^ - -j  =0  v  e  E(Q). 

The  principle  of  virtual  work  in  the  case  of  initial  stress: 

/  CijkiUijVk,i  dV  +  /  cr9.u  v  •  dV  = 

Ja  Ja  J 

f  FiVi  dV  +  f  T{Vi  dS  +  /  T CijkiakiVij  dV 
J  3£lrr 
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Ld.Fct=  1.8346e-01 
Run  =  4 
Fnc.  =  Formula 
Max  =  4.5846e-02 
Min  =  -4.5847e-02 

4.5000e-02 

3.5000e-02 

2.5000e-02 

1.5000e-02 

5.0000e-03 

-5.0000e-03 

-1.5000e-02 

-2.5000e-02 

-3.5000e-02 

-4.5000e-02 


Figure  12:  Lockheed  test  problem  2:  The  function  un  =  uxnx  +  uyny 


for  all  Vi  £  E(tt). 

Although  two  mathematical  models  of  the  general  theory  of  elastic  stability  exist  in  the 
classical  literature,  some  fundamental  questions  concerning  the  existence  of  a  solution,  the 
properties  of  the  spectrum,  and  their  relationship  to  loss  of  stability  had  not  been  investigated 
previously.  Two  working  hypotheses  were  advanced: 

1.  The  spectrum  is  a  point  spectrum,  hence  it  is  meaningful  to  consider  the  lowest  nonzero 
eigenvalue  as  an  indicator  of  the  onset  of  instability; 

2.  The  minimal  eigenvalue  of  the  finite  dimensional  problem  converges  to  its  infinite 
dimensional  counterpart  as  the  finite  element  space  is  enlarged  (i.e.,  the  degrees  of 
freedom  are  increased). 
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Both  classical  formulations  were  implemented  in  fully  three-dimensional  setting  so  that  nu¬ 
merical  experiments  could  performed.  For  thin  structures  the  results  closely  matched  the 
classical  results.  It  was  found  that  the  two  classical  models  yield  virtually  identical  results. 

In  a  parallel  investigation  the  first  working  hypothesis  was  proven  by  Professors  Manil 
Suri  and  Ivo  Babuska.  At  present  it  is  not  known  whether  the  second  hypothesis  can  be 
proven,  but  the  available  numerical  results  have  not  contradicted  it. 

The  relationship  between  the  limits  of  elastic  stability  estimated  by  the  use  of  linear 
models  and  incremental  models  was  investigated.  It  was  found  that  for  conservative  loads 
a  close  relationship  exists  but  the  treatment  of  follower  loads  through  the  solution  of  linear 
eigenvalue  problems  does  not  appear  possible,  with  the  exception  of  very  special  cases,  such 
as  the  buckling  of  circular  rings.  See,  for  example,  [9]. 

The  problem  of  modeling  the  elastic  buckling  of  fibers  in  fiber-matrix  composites  was 
investigated.  A  model  problem  has  been  solved.  For  the  investigation  of  the  stability  of  a 
large  number  of  fibers  the  use  of  periodic  boundary  conditions  is  necessary.  Implementation 
of  periodic  boundary  conditions  and  further  investigation  of  the  stability  of  fibers  is  being 
planned  for  1998. 
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